#include <kuai/misc/thread.hpp>
#include <kuai/misc/filter.hpp>
#include <kuai/sim/energy.hpp>
#include <kuai/sim/efunctor.hpp>


#ifndef _KUAI_SIM_NONBOND_
#define _KUAI_SIM_NONBOND_

namespace kuai {

	
	class BasicNonbondEnergy
		: public BasicForceFieldEnergy
	{
	public:
		enum {
			MAX_ATOMS_IN_CHARGE_GROUP = 64
		};

		BasicNonbondEnergy(Index rank0)
			: BasicForceFieldEnergy(rank0)
		{ 
			k_columnb = NULL;
			charges = NULL;
			atomType = NULL;
			groupID = NULL;
			cutoff = 0;
			nAtoms = 0;
			par0 = 0;
		}
		RealNumber* k_columnb;
		RealNumber* charges;
		Index* atomType;
		Index* groupID;
		Byte* filter;
		Index par0;
		RealNumber cutoff;
		Index nAtoms;
		RealNumber buffer;

		bool test(size_t a1, size_t a2) const {
			static const size_t BITS_OF_FILTER = 8 * sizeof(Byte);
			static const Char mask[] = {1, 2, 4, 8, 16, 32, 64, char(128)};
			assert (a1 > a2);
			size_t offset = a1 * (a1-1) / 2 + a2;
			size_t m = offset / BITS_OF_FILTER, n = offset % BITS_OF_FILTER;
			return (filter[m] & mask[n]) != 0;
		}
	};


	template <typename Functor, bool addColumnb, bool usePBC>
		class NonbondEnergy
			: public BasicNonbondEnergy
	{
	public: // Interface to EnergyFunction
		NonbondEnergy()
			: BasicNonbondEnergy(addColumnb?2:1)
		{ 
			pFuncAddE = NULL;
			pFuncAddEG = NULL;
			pFuncAddEGH = NULL;
		}

		virtual void start(SimulationEnvironment& env) {
			if (env.isMainThread()) {
				if (k_columnb) {
					columnb_const = *k_columnb;
				}
				if (nlist > 0) {
					pFuncAddE = &NonbondEnergy::addE_List;
					pFuncAddEG = &NonbondEnergy::addEG_List;
					pFuncAddEGH = &NonbondEnergy::addEGH_List;
					pFuncAddHX = &NonbondEnergy::addHX_List;
				}
				else if (cutoff <= 0) {
					pFuncAddE = &NonbondEnergy::addE_NoC;
					pFuncAddEG = &NonbondEnergy::addEG_NoC;
					pFuncAddEGH = &NonbondEnergy::addEGH_NoC;
					pFuncAddHX = &NonbondEnergy::addHX_NoC;
				}
				else if (groupID == NULL) {
					pFuncAddE = &NonbondEnergy::addE_ABC;
					pFuncAddEG = &NonbondEnergy::addEG_ABC;
					pFuncAddEGH = &NonbondEnergy::addEGH_ABC;
					pFuncAddHX = &NonbondEnergy::addHX_ABC;
				}
				else {
					pFuncAddE = &NonbondEnergy::addE_GBC;
					pFuncAddEG = &NonbondEnergy::addEG_GBC;
					pFuncAddEGH = &NonbondEnergy::addEGH_GBC;
					pFuncAddHX = &NonbondEnergy::addHX_GBC;
				}
			}
			env.wait();
			if (pFuncAddE == &NonbondEnergy::addE_ABC) {
				if (env.isMainThread()) {
					coords0.resize(nAtoms);	 // prepare the coords0 
					bufferList.reserve(64*1024*3);	// 
				}
				env.wait();

				updateABCBuffer(env);
			}
			else if (pFuncAddE == &NonbondEnergy::addE_GBC) {
				prepareGBC(env);
			}
		}

		virtual void addE(SimulationEnvironment& env, RealNumber e[]) {
			(this->*pFuncAddE)(env, e);
		}
		virtual void addEG(SimulationEnvironment& env, RealNumber e[], XYZ gred[]) {
			(this->*pFuncAddEG)(env, e, gred);
		}
		virtual void addEGH(SimulationEnvironment& env, RealNumber e[], XYZ gred[], RealNumber* hess[]) {
			(this->*pFuncAddEGH)(env, e, gred, hess);
		}
		virtual void addHX(SimulationEnvironment& env, const RealNumber x[], RealNumber result[]) {
			(this->*pFuncAddHX)(env, x, result);
		}

	
	private:
		void addE_List(SimulationEnvironment& env, Index* lists, size_t nlist, RealNumber e[]) {
			assert (nlist % 3 == 0);
			SizePair range = env.range(nlist / 3);
			for (size_t i = range.first; i != range.second; ++i) {
				Index atom1 = lists[i*3];
				Index atom2 = lists[i*3+1];
				Index iPar  = lists[i*3+2];
				addEItem(atom1, atom2, iPar,e);
			}
		}
		void addE_List(SimulationEnvironment& env, RealNumber e[]) {
			addE_List(env, lists, nlist, e);
		}

		void addEG_List(SimulationEnvironment& env, Index* lists, size_t nlist, RealNumber e[], XYZ gred[]) {
			assert (nlist % 3 == 0);
			SizePair range = env.range(nlist / 3);
			for (size_t i = range.first; i != range.second; ++i) {
				Index atom1 = lists[i*3];
				Index atom2 = lists[i*3+1];
				Index iPar  = lists[i*3+2];
				addEGItem(atom1, atom2, iPar, e, gred);
			}
		}
		void addEGH_List(SimulationEnvironment& env, Index* lists, size_t nlist, RealNumber e[], XYZ gred[], RealNumber* hess[]) {
			assert (nlist % 3 == 0);
			SizePair range = env.range(nlist / 3);
			for (size_t i = range.first; i != range.second; ++i) {
				Index atom1 = lists[i*3];
				Index atom2 = lists[i*3+1];
				Index iPar  = lists[i*3+2];
				addEGHItem(env, atom1, atom2, iPar, e, gred, hess);
			}
		}
		void addHX_List(SimulationEnvironment& env, Index* lists, size_t nlist, const RealNumber x[], RealNumber result[]) {
			assert (nlist % 3 == 0);
			SizePair range = env.range(nlist / 3);
			for (size_t i = range.first; i != range.second; ++i) {
				Index atom1 = lists[i*3];
				Index atom2 = lists[i*3+1];
				Index iPar  = lists[i*3+2];
				addHXItem(atom1, atom2, iPar, x, result);
			}
		}

		void addEG_List(SimulationEnvironment& env, RealNumber e[], XYZ gred[]) {
			addEG_List(env, lists, nlist, e, gred);
		}
		void addEGH_List(SimulationEnvironment& env, RealNumber e[], XYZ gred[], RealNumber* hess[]) {
			addEGH_List(env, lists, nlist, e, gred, hess);
		}
		void addHX_List(SimulationEnvironment& env, const RealNumber x[], RealNumber result[]) {
			addHX_List(env, lists, nlist, x, result);
		}

		void addE_NoC(SimulationEnvironment& env, RealNumber e[]) {
			SizePair range = env.range(nAtoms/2);
			for (size_t atom1 = range.first; atom1 != range.second; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEItem(atom1, atom2, iPar,e);
					}
				}
			}
			if (nAtoms % 2 == 1) {
				size_t atom1 = (nAtoms / 2) + 1;
				std::pair<size_t, size_t> range = env.range(atom1);
				for (size_t atom2 = range.first; atom2 < range.second; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEItem(atom1, atom2, iPar,e);
					}
				}
			}
			for (size_t atom1 = nAtoms - range.second; atom1 != nAtoms - range.first; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEItem(atom1, atom2, iPar,e);
					}
				}
			}
		}

		void addEG_NoC(SimulationEnvironment& env, RealNumber e[], XYZ gred[]) {
			SizePair range = env.range(nAtoms/2);
			for (size_t atom1 = range.first; atom1 != range.second; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEGItem(atom1, atom2, iPar, e, gred);
					}
				}
			}
			if (nAtoms % 2 == 1) {
				size_t atom1 = (nAtoms / 2) + 1;
				std::pair<size_t, size_t> range = env.range(atom1);
				for (size_t atom2 = range.first; atom2 < range.second; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEGItem(atom1, atom2, iPar, e, gred);
					}
				}
			}
			for (size_t atom1 = nAtoms - range.second; atom1 != nAtoms - range.first; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEGItem(atom1, atom2, iPar, e, gred);
					}
				}
			}
		}

		void addEGH_NoC(SimulationEnvironment& env, RealNumber e[], XYZ gred[], RealNumber* hess[]) {
			SizePair range = env.range(nAtoms/2);
			for (size_t atom1 = range.first; atom1 != range.second; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEGHItem(env, atom1, atom2, iPar, e, gred, hess);
					}
				}
			}
			if (nAtoms % 2 == 1) {
				size_t atom1 = (nAtoms / 2) + 1;
				std::pair<size_t, size_t> range = env.range(atom1);
				for (size_t atom2 = range.first; atom2 < range.second; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEGHItem(env, atom1, atom2, iPar, e, gred, hess);
					}
				}
			}
			for (size_t atom1 = nAtoms - range.second; atom1 != nAtoms - range.first; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addEGHItem(env, atom1, atom2, iPar, e, gred, hess);
					}
				}
			}
		}
		void addHX_NoC(SimulationEnvironment& env, const RealNumber x[], RealNumber result[]) {
			SizePair range = env.range(nAtoms/2);
			for (size_t atom1 = range.first; atom1 != range.second; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addHXItem(atom1, atom2, iPar, x, result);
					}
				}
			}
			if (nAtoms % 2 == 1) {
				size_t atom1 = (nAtoms / 2) + 1;
				std::pair<size_t, size_t> range = env.range(atom1);
				for (size_t atom2 = range.first; atom2 < range.second; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addHXItem(atom1, atom2, iPar, x, result);
					}
				}
			}
			for (size_t atom1 = nAtoms - range.second; atom1 != nAtoms - range.first; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					if (test(atom1, atom2)) {
						Index iPar  = getParIndex(atom1, atom2);
						addHXItem(atom1, atom2, iPar, x, result);
					}
				}
			}
		}

		void addE_ABC(SimulationEnvironment& env, RealNumber e[]) {
			if (needUpdateABCBuffer(env)) {
				updateABCBuffer(env);
			}
			std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
			RealNumber cutoff2=cutoff*cutoff;
			for (size_t i = range.first; i < range.second; ++i) {
				Index atom1 = bufferList[i*3];
				Index atom2 = bufferList[i*3+1];
				Index iPar  = bufferList[i*3+2];
				XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
				if (usePBC) {
					pbc->update(delta);
				}
				RealNumber r2;
				if (shorter2(delta, cutoff2, r2)) {
					addEItem(sqrt(r2), atom1, atom2, iPar,e);
				}
			}
		}

		void addEG_ABC(SimulationEnvironment& env, RealNumber e[], XYZ gred[]) {
			if (needUpdateABCBuffer(env)) {
				updateABCBuffer(env);
			}
			std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
			RealNumber cutoff2=cutoff*cutoff;
			for (size_t i = range.first; i < range.second; ++i) {
				Index atom1 = bufferList[i*3];
				Index atom2 = bufferList[i*3+1];
				Index iPar  = bufferList[i*3+2];
				XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
				if (usePBC) {
					pbc->update(delta);
				}
				RealNumber r2;
				if (shorter2(delta, cutoff2, r2)) {
					addEGItem(delta, sqrt(r2), atom1, atom2, iPar, e, gred);
				}
			}
		}

		void addEGH_ABC(SimulationEnvironment& env, RealNumber e[], XYZ gred[], RealNumber* hess[]) {
			if (needUpdateABCBuffer(env)) {
				updateABCBuffer(env);
			}
			std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
			RealNumber cutoff2=cutoff*cutoff;
			for (size_t i = range.first; i < range.second; ++i) {
				Index atom1 = bufferList[i*3];
				Index atom2 = bufferList[i*3+1];
				Index iPar  = bufferList[i*3+2];
				XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
				if (usePBC) {
					pbc->update(delta);
				}
				RealNumber r2;
				if (shorter2(delta, cutoff2, r2)) {
					addEGHItem(env, delta, atom1, atom2, iPar, e, gred, hess);
				}
			}
		}
		void addHX_ABC(SimulationEnvironment& env, const RealNumber x[], RealNumber result[]) {
			if (needUpdateABCBuffer(env)) {
				updateABCBuffer(env);
			}
			std::pair<size_t, size_t> range = env.range(bufferList.size()/3);
			RealNumber cutoff2=cutoff*cutoff;
			for (size_t i = range.first; i < range.second; ++i) {
				Index atom1 = bufferList[i*3];
				Index atom2 = bufferList[i*3+1];
				Index iPar  = bufferList[i*3+2];
				XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
				if (usePBC) {
					pbc->update(delta);
				}
				RealNumber r2;
				if (shorter2(delta, cutoff2, r2)) {
					addHXItem(delta, atom1, atom2, iPar, x, result);
				}
			}
		}

		void addE_GBC(SimulationEnvironment& env, RealNumber e[]) {
			if (needUpdateGBCBuffer(env)) {
				updateGBCBuffer(env);
			}
			addE_List(env, &innerGroupList[0], innerGroupList.size(), e);

			if (env.isMainThread()) {
				count = 0;
			}
			env.wait();
			Index nGroup = groupCenters.size();
			Index i = nGroup - nextCount(env) - 1;
			RealNumber r2 = cutoff * cutoff;
			while (0 < i && i < nGroup) {
				const XYZ& coordI = coords[groupCenters[i]];
				for (size_t j = 0; j < i; ++j) {
					const XYZ& coordJ = coords[groupCenters[j]];
					XYZ delta(coordI.x - coordJ.x, coordI.y - coordJ.y, coordI.z - coordJ.z);
					if (usePBC) {
						pbc->update(delta);
					}
					if (shorter2(delta, r2)) {
						addE_GBC_Item(i, j, e);
					}
				}
				i = nGroup - nextCount(env) - 1;
			}
		}

		void addEG_GBC(SimulationEnvironment& env, RealNumber e[], XYZ gred[]) {
			if (needUpdateGBCBuffer(env)) {
				updateGBCBuffer(env);
			}
			addEG_List(env, &innerGroupList[0], innerGroupList.size(), e, gred);

			if (env.isMainThread()) {
				count = 0;
			}
			env.wait();
			Index nGroup = groupCenters.size();
			Index i = nGroup - nextCount(env) - 1;
			RealNumber r2 = cutoff * cutoff;
			while (0 < i && i < nGroup) {
				const XYZ& coordI = coords[groupCenters[i]];
				for (size_t j = 0; j < i; ++j) {
					const XYZ& coordJ = coords[groupCenters[j]];
					XYZ delta(coordI.x - coordJ.x, coordI.y - coordJ.y, coordI.z - coordJ.z);
					if (usePBC) {
						pbc->update(delta);
					}
					if (shorter2(delta, r2)) {
						addEG_GBC_Item(i, j, e, gred);
					}
				}
				i = nGroup - nextCount(env) - 1;
			}
		}

		void addEGH_GBC(SimulationEnvironment& env, RealNumber e[], XYZ gred[], RealNumber* hess[]) {
			if (needUpdateGBCBuffer(env)) {
				updateGBCBuffer(env);
			}
			addEGH_List(env, &innerGroupList[0], innerGroupList.size(), e, gred, hess);

			if (env.isMainThread()) {
				count = 0;
			}
			env.wait();
			Index nGroup = groupCenters.size();
			Index i = nGroup - nextCount(env) - 1;
			RealNumber r2 = cutoff * cutoff;
			while (0 < i && i < nGroup) {
				const XYZ& coordI = coords[groupCenters[i]];
				for (size_t j = 0; j < i; ++j) {
					const XYZ& coordJ = coords[groupCenters[j]];
					XYZ delta(coordI.x - coordJ.x, coordI.y - coordJ.y, coordI.z - coordJ.z);
					if (usePBC) {
						pbc->update(delta);
					}
					if (shorter2(delta, r2)) {
						addEGH_GBC_Item(env, i, j, e, gred, hess);
					}
				}
				i = nGroup - nextCount(env) - 1;
			}
		}

		void addHX_GBC(SimulationEnvironment& env, const RealNumber x[], RealNumber result[]) {
			if (needUpdateGBCBuffer(env)) {
				updateGBCBuffer(env);
			}
			addHX_List(env, &innerGroupList[0], innerGroupList.size(), x, result);

			if (env.isMainThread()) {
				count = 0;
			}
			env.wait();
			Index nGroup = groupCenters.size();
			Index i = nGroup - nextCount(env) - 1;
			RealNumber r2 = cutoff * cutoff;
			while (0 < i && i < nGroup) {
				const XYZ& coordI = coords[groupCenters[i]];
				for (size_t j = 0; j < i; ++j) {
					const XYZ& coordJ = coords[groupCenters[j]];
					XYZ delta(coordI.x - coordJ.x, coordI.y - coordJ.y, coordI.z - coordJ.z);
					if (usePBC) {
						pbc->update(delta);
					}
					if (shorter2(delta, r2)) {
						addHX_GBC_Item(env, i, j, x, result);
					}
				}
				i = nGroup - nextCount(env) - 1;
			}
		}


		void addE_GBC_Item(Index i, Index j, RealNumber e[]) const {
			for (size_t ii = 0; ii < groupAtomCount[i]; ++ii) {
				for (size_t jj = 0; jj < groupAtomCount[j]; ++jj) {
					Index atomI = groupAtoms[i*MAX_ATOMS_IN_CHARGE_GROUP+ii];
					Index atomJ = groupAtoms[j*MAX_ATOMS_IN_CHARGE_GROUP+jj];
					if (atomI < atomJ) {
						std::swap(atomI, atomJ);
					}
					if (test(atomI, atomJ)) {
						Index iPar = getParIndex(atomI, atomJ);
						addEItem(atomI, atomJ, iPar, e);
					}
				}
			}
		}

		void addEG_GBC_Item(Index i, Index j, RealNumber e[], XYZ gred[]) const {
			for (size_t ii = 0; ii < groupAtomCount[i]; ++ii) {
				for (size_t jj = 0; jj < groupAtomCount[j]; ++jj) {
					Index atomI = groupAtoms[i*MAX_ATOMS_IN_CHARGE_GROUP+ii];
					Index atomJ = groupAtoms[j*MAX_ATOMS_IN_CHARGE_GROUP+jj];
					if (atomI < atomJ) {
						std::swap(atomI, atomJ);
					}
					if (test(atomI, atomJ)) {
						Index iPar = getParIndex(atomI, atomJ);
						addEGItem(atomI, atomJ, iPar, e, gred);
					}
				}
			}
		}

		void addEGH_GBC_Item(SimulationEnvironment& env, Index i, Index j, RealNumber e[], XYZ gred[], RealNumber* hess[]) const {
			for (size_t ii = 0; ii < groupAtomCount[i]; ++ii) {
				for (size_t jj = 0; jj < groupAtomCount[j]; ++jj) {
					Index atomI = groupAtoms[i*MAX_ATOMS_IN_CHARGE_GROUP+ii];
					Index atomJ = groupAtoms[j*MAX_ATOMS_IN_CHARGE_GROUP+jj];
					if (atomI < atomJ) {
						std::swap(atomI, atomJ);
					}
					if (test(atomI, atomJ)) {
						Index iPar = getParIndex(atomI, atomJ);
						addEGHItem(env, atomI, atomJ, iPar, e, gred, hess);
					}
				}
			}
		}

		void addHX_GBC_Item(SimulationEnvironment& env, Index i, Index j, const RealNumber x[], RealNumber result[]) const {
			for (size_t ii = 0; ii < groupAtomCount[i]; ++ii) {
				for (size_t jj = 0; jj < groupAtomCount[j]; ++jj) {
					Index atomI = groupAtoms[i*MAX_ATOMS_IN_CHARGE_GROUP+ii];
					Index atomJ = groupAtoms[j*MAX_ATOMS_IN_CHARGE_GROUP+jj];
					if (atomI < atomJ) {
						std::swap(atomI, atomJ);
					}
					if (test(atomI, atomJ)) {
						Index iPar = getParIndex(atomI, atomJ);
						addHXItem(atomI, atomJ, iPar, x, result);
					}
				}
			}
		}

		Index getParIndex(size_t a1, size_t a2) const {
			Index t1 = atomType[a1], t2 = atomType[a2];
			if (t1 < t2) {
				std::swap(t1, t2);
			}
			Index result = par0 + (t1 * (t1+1) / 2 + t2) * Functor::PAR_NUMBER;
			return result;
		}
		void addEItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[]) const {
			addEItem(abs(delta), i1, i2, iPar, e);
		}
		void addEItem(RealNumber r, Index i1, Index i2, Index iPar, RealNumber e[]) const {
			RealNumber rr = 1/r;
			e[0] += func(r, pars+iPar);

			if (addColumnb) {
				e[1] += columnb_const * charges[i1]*charges[i2] * rr;
			}
		}
		void addEItem(Index i1, Index i2, Index iPar, RealNumber e[]) const {
			XYZ delta(coords[i1].x - coords[i2].x, coords[i1].y - coords[i2].y, coords[i1].z - coords[i2].z);
			if (usePBC) {
				pbc->update(delta);
			}
			addEItem(delta, i1, i2, iPar, e);
		}

		void addEGItem(Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const
		{
			XYZ delta(coords[i1].x - coords[i2].x, coords[i1].y - coords[i2].y, coords[i1].z - coords[i2].z);
			if (usePBC) {
				pbc->update(delta);
			}
			addEGItem(delta, i1, i2, iPar, e, gred);
		}
		void addEGHItem(SimulationEnvironment& env, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[], RealNumber* hess[]) const
		{
			XYZ delta(coords[i1].x - coords[i2].x, coords[i1].y - coords[i2].y, coords[i1].z - coords[i2].z);
			if (usePBC) {
				pbc->update(delta);
			}
			addEGHItem(env, delta, i1, i2, iPar, e, gred, hess);
		}
		void addHXItem(Index i1, Index i2, Index iPar, const RealNumber x[], RealNumber result[]) const
		{
			XYZ delta(coords[i1].x - coords[i2].x, coords[i1].y - coords[i2].y, coords[i1].z - coords[i2].z);
			if (usePBC) {
				pbc->update(delta);
			}
			addHXItem(delta, i1, i2, iPar, x, result);
		}

		void addEGItem(XYZ& delta, RealNumber r, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const
		{
			RealNumber rr = 1/r;
			RealNumber dedr;
			e[0] += func(r, pars+iPar, dedr);

			if (addColumnb) {
				RealNumber k = columnb_const * charges[i1]*charges[i2] * rr;
				e[1] += k;
				dedr -= k * rr;
			}
			delta *= dedr * rr;
			gred[i1] += delta;
			gred[i2] -= delta;
		}

		void addEGItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const
		{
			addEGItem(delta, abs(delta), i1, i2, iPar, e, gred);
		}
		void addEGHItem(SimulationEnvironment& env, const XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[], RealNumber* hess[]) const
		{
			XYZ dx[1];
			HessianItem dx2[2];

			RealNumber r = calcAtomPair(delta, dx, dx2);
			RealNumber rr = 1/r;
			RealNumber dedr, d2edr2;
			e[0] += func(r, pars+iPar, dedr, d2edr2);

			if (addColumnb) {
				RealNumber k = columnb_const * charges[i1]*charges[i2] * rr;
				e[1] += k;
				k *= rr;
				dedr -= k;
				k *= 2 * rr;
				d2edr2 += k;
			}

			dx2[0] *= dedr;
			dx2[1] = prod(dx[0], dx[0]);
			dx2[1] *= d2edr2;
			dx2[0] += dx2[1];

			dx[0] *= dedr;
			gred[i1] += dx[0];
			gred[i2] -= dx[0];

			env.lock();
				addHessianItem(dx2[0], i1, i1, hess);
				subHessianItem(dx2[0], i1, i2, hess);
				subHessianItem(dx2[0], i2, i1, hess);
				addHessianItem(dx2[0], i2, i2, hess);
			env.free();
		}

		void addHXItem(const XYZ& delta, Index i1, Index i2, Index iPar, const RealNumber x[], RealNumber result[]) const
		{
			XYZ dx[1];
			HessianItem dx2[2];

			RealNumber r = calcAtomPair(delta, dx, dx2);
			RealNumber rr = 1/r;
			RealNumber dedr, d2edr2;
			func(r, pars+iPar, dedr, d2edr2);

			if (addColumnb) {
				RealNumber k = columnb_const * charges[i1]*charges[i2] * rr;
				k *= rr;
				dedr -= k;
				k *= 2 * rr;
				d2edr2 += k;
			}

			dx2[0] *= dedr;
			dx2[1] = prod(dx[0], dx[0]);
			dx2[1] *= d2edr2;
			dx2[0] += dx2[1];

			addHessianItem(dx2[0], i1, i1, x, result);
			subHessianItem(dx2[0], i1, i2, x, result);
			subHessianItem(dx2[0], i2, i1, x, result);
			addHessianItem(dx2[0], i2, i2, x, result);
		}

		bool needUpdateABCBuffer(SimulationEnvironment& env) {
			std::pair<size_t, size_t> range = env.range(nAtoms);
			RealNumber r2 = buffer*buffer / 4; // square of half buffer

			for (size_t i = range.first; !needUpdate && i < range.second; ++i) {
				XYZ delta(coords0[i].x-coords[i].x, coords0[i].y-coords[i].y, coords0[i].z-coords[i].z);
				if (!shorter2(delta, r2)) {
					needUpdate = true;
					break;
				}
			}
			env.wait();
			return needUpdate;
		}

		void updateABCBuffer(SimulationEnvironment& env, Index atom1, Index atom2, RealNumber r2) {
			if (test(atom1, atom2)) {
				XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
				if (usePBC) {
					pbc->update(delta);
				}
				if (shorter2(delta, r2)) {
					env.lock();
					bufferList.push_back(atom1);
					bufferList.push_back(atom2);
					bufferList.push_back(getParIndex(atom1, atom2));
					env.free();
				}
			}
		}

		void updateABCBuffer(SimulationEnvironment& env) {
			if (env.isMainThread()) {
				bufferList.clear();
			}
			env.wait();
			std::pair<size_t, size_t> range = env.range(nAtoms);
			std::memcpy(&coords0[range.first], &coords[range.first], sizeof(XYZ)*(range.second-range.first));

			RealNumber r2 = (cutoff+buffer)*(cutoff+buffer);
			range = env.range(nAtoms/2);
			for (size_t atom1 = range.first; atom1 < range.second; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					updateABCBuffer(env, atom1, atom2, r2);
				}
			}
			if (nAtoms % 2 == 1) {
				size_t atom1 = (nAtoms / 2) + 1;
				std::pair<size_t, size_t> range = env.range(atom1);
				for (size_t atom2 = range.first; atom2 < range.second; ++atom2) {
					updateABCBuffer(env, atom1, atom2, r2);
				}
			}
			for (size_t atom1 = nAtoms - range.second; atom1 != nAtoms - range.first; ++atom1) {
				for (size_t atom2 = 0; atom2 < atom1; ++atom2) {
					updateABCBuffer(env, atom1, atom2, r2);
				}
			}
			if (env.isMainThread()) {
				needUpdate = false;
			}
			env.wait();
		}

		
		void prepareGBC(SimulationEnvironment& env) {
			makeUpGroup(env);
			makeUpInnerGroupList(env);
			if (env.isMainThread()) {
				coords0.resize(groupCenters.size());
			} 
			env.wait();
			updateGBCBuffer(env);
		}

		void makeUpGroup(SimulationEnvironment& env) {
			if (env.isMainThread()) {
				count = 0;
				groupCenters.clear();
				groupAtomCount.clear();
				groupAtoms.clear();
			}
			env.wait();
			std::pair<size_t, size_t> range = env.range(nAtoms);
			Index nGroup = 0;
			for (size_t i = range.first; i < range.second; ++i) {
				if (groupID[i] > nGroup) {
					nGroup = groupID[i];
				}
			}
			nGroup += 1;
			env.lock();
			if (groupCenters.size() < nGroup) {
				groupCenters.resize(nGroup);
				groupAtomCount.resize(nGroup, 0);
				groupAtoms.resize(nGroup*MAX_ATOMS_IN_CHARGE_GROUP);
			}
			env.free();
			env.wait();

			for (size_t i = range.first; i < range.second; ++i) {
				Index gid = groupID[i];
				env.lock();
				Index offset = groupAtomCount[gid]++;
				env.free();
				assert (offset < MAX_ATOMS_IN_CHARGE_GROUP);
				offset += MAX_ATOMS_IN_CHARGE_GROUP*gid;
				groupAtoms[offset] = i;
			}
			env.wait();

			nGroup = groupCenters.size();

			Index gid = nextCount(env);
			while (gid < nGroup) {
				Index n = groupAtomCount[gid];

				if (n == 1 || charges == NULL) {
					groupCenters[gid] = groupAtoms[gid*MAX_ATOMS_IN_CHARGE_GROUP];
				}
				else if (n > 1) {
					XYZ center(0, 0, 0);
					RealNumber totalCharge = 0;
					for (size_t j = 0; j < n; ++j) {
						Index atomJ = groupAtoms[gid*MAX_ATOMS_IN_CHARGE_GROUP+j];
						totalCharge += fabs(charges[atomJ]);
						center += fabs(charges[atomJ])*coords[atomJ];
					}
					Index atomI = groupAtoms[gid*MAX_ATOMS_IN_CHARGE_GROUP];
					groupCenters[gid] = atomI;
					if (totalCharge > 1e-5) {
						center /= totalCharge;
						XYZ delta = coords[atomI] - center;
						if (usePBC) {
							pbc->update(delta);
						}
						RealNumber rI2 = dot(delta, delta);
						for (size_t j = 1; j < n; ++j) {
							Index atomJ = groupAtoms[gid*MAX_ATOMS_IN_CHARGE_GROUP+j];
							delta = coords[atomJ] - center;
							if (usePBC) {
								pbc->update(delta);
							}
							RealNumber rJ2 = dot(delta, delta);
							if (rJ2 < rI2) {
								rI2 = rJ2;
								groupCenters[gid] = atomJ;
							}
						}
					}
				}
				else {
					assert (n == 0);
				}
				
				gid = nextCount(env);
			}
			env.wait();
		}

		void makeUpInnerGroupList(SimulationEnvironment& env) {
			if (env.isMainThread()) {
				count = 0;
				innerGroupList.clear();
			}
			Index nGroup = groupCenters.size();
			Index gid = nextCount(env);
			while (gid < nGroup) {
				Index n = groupAtomCount[gid];
				const Index* pStart = &groupAtoms[gid*MAX_ATOMS_IN_CHARGE_GROUP];

				if (n > 1) {
					for (size_t i = 1; i < n; ++i) {
						for (size_t j = 0; j < i; ++j) {
							Index atomI = pStart[i];
							Index atomJ = pStart[j];
							if (atomI < atomJ) {
								std::swap(atomI, atomJ);
							}
							if (test(atomI, atomJ)) {
								Index iPar = getParIndex(atomI, atomJ);
								env.lock();
								innerGroupList.push_back(atomI);
								innerGroupList.push_back(atomJ);
								innerGroupList.push_back(iPar);
								env.free();
							}
						}
					}
				}
				
				gid = nextCount(env);
			}
		
		}


		bool needUpdateGBCBuffer(SimulationEnvironment& env) {
			std::pair<size_t, size_t> range = env.range(groupCenters.size());
			RealNumber r2 = buffer*buffer / 4; // square of half buffer

			for (size_t i = range.first; !needUpdate && i < range.second; ++i) {
				Index atomI = groupCenters[i];
				XYZ delta(coords0[i].x-coords[atomI].x, coords0[i].y-coords[atomI].y, coords0[i].z-coords[atomI].z);
				if (!shorter2(delta, r2)) {
					needUpdate = true;
					break;
				}
			}
			env.wait();
			return needUpdate;
		}
		void updateGBCBuffer(SimulationEnvironment& env, Index gi, Index gj, RealNumber r2) {
			Index atom1 = groupCenters[gi], atom2 = groupCenters[gj];
			XYZ delta(coords[atom1].x - coords[atom2].x, coords[atom1].y - coords[atom2].y, coords[atom1].z - coords[atom2].z);
			if (usePBC) {
				pbc->update(delta);
			}
			if (shorter2(delta, r2)) {
				env.lock();
				bufferList.push_back(gi);
				bufferList.push_back(gj);
				env.free();
			}
		}
		void updateGBCBuffer(SimulationEnvironment& env) {
			if (env.isMainThread()) {
				bufferList.clear();
			}
			env.wait();
			
			Index nGroups = groupCenters.size();
			RealNumber r2 = (cutoff+buffer)*(cutoff+buffer);
			std::pair<size_t, size_t> range = env.range(nGroups/2);
			for (size_t i = range.first; i < range.second; ++i) {
				coords0[i] = coords[groupCenters[i]];
				for (size_t j = 0; j < i; ++j) {
					updateGBCBuffer(env, i, j, r2);
				}
			}
			if (nGroups % 2 == 1) {
				size_t i = (nGroups / 2) + 1;
				coords0[i] = coords[groupCenters[i]];
				std::pair<size_t, size_t> range = env.range(i);
				for (size_t j = range.first; j < range.second; ++j) {
					updateGBCBuffer(env, i, j, r2);
				}
			}
			for (size_t i = nGroups - range.second; i != nGroups - range.first; ++i) {
				coords0[i] = coords[groupCenters[i]];
				for (size_t j = 0; j < i; ++j) {
					updateGBCBuffer(env, i, j, r2);
				}
			}
			
			if (env.isMainThread()) {
				needUpdate = false;
			}
			env.wait();
		}

		Index nextCount(SimulationEnvironment& env) {
			env.lock();
			Index result = count++;
			env.free();
			return result;
		}
		
	private:
		bool needUpdate;
		std::vector<Index> bufferList;
		std::vector<Index> innerGroupList;
		std::vector<XYZ> coords0;

		std::vector<Index> groupCenters;
		std::vector<Index> groupAtomCount;
		std::vector<Index> groupAtoms;

		Index count;

		RealNumber columnb_const;
		Functor func;
		void (NonbondEnergy::*pFuncAddE)(SimulationEnvironment&, RealNumber []);
		void (NonbondEnergy::*pFuncAddEG)(SimulationEnvironment&, RealNumber [], XYZ []);
		void (NonbondEnergy::*pFuncAddHX)(SimulationEnvironment&, const RealNumber x[], RealNumber result[]);
		void (NonbondEnergy::*pFuncAddEGH)(SimulationEnvironment&, RealNumber [], XYZ [], RealNumber* hess[]);
	};

	BasicNonbondEnergy* createNonbondEnergyFunction(const char szName[], 
		const PBC* pbc, bool hasCharge);

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[]) const;
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, false>::addEItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[]) const;

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEGItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const;
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, false>::addEGItem(XYZ& delta, Index i1, Index i2, Index iPar, RealNumber e[], XYZ gred[]) const;

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addE_ABC(SimulationEnvironment& env, RealNumber e[]);
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, false>::addE_ABC(SimulationEnvironment& env, RealNumber e[]);

	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEG_ABC(SimulationEnvironment& env, RealNumber e[], XYZ gred[]);
	template <>
		void NonbondEnergy<FunctorLJ12_6, false, true>::addEG_ABC(SimulationEnvironment& env, RealNumber e[], XYZ gred[]);
	
	
	
}

#endif
